Recursive oo Group lasso 

Yilun Chen, Student Member, IEEE, and Alfred O. Hero, III, Fellow, IEEE 



Abstract — We introduce a recursive adaptive group lasso al- 
gorithm for real-time penalized least squares prediction that 
produces a time sequence of optimal sparse predictor coefficient 
vectors. At each time index the proposed algorithm computes an 
exact update of the optimal £1,00 -penalized recursive least squares 
(RLS) predictor. Each update minimizes a convex but non- 
differentiable function optimization problem. We develop an on- 
line homotopy method to reduce the computational complexity. 
Numerical simulations demonstrate that the proposed algorithm 
outperforms the li regularized RLS algorithm for a group sparse 
system identification problem and has lower implementation 
complexity than direct group lasso solvers. 

Index Terms — RLS, group sparsity, mixed norm, homotopy, 
group lasso, system identification 

I. Introduction 

Recursive Least Squares (RLS) is a widely used method for 
adaptive filtering and prediction in signal processing and re- 
lated fields. Its applications include: acoustic echo cancelation; 
wireless channel equalization; interference cancelation and 
data streaming predictors. In these applications a measurement 
stream is recursively fitted to a linear model, described by 
the coefficients of an FIR prediction filter, in such a way to 
minimize a weighted average of squared residual prediction 
errors. Compared to other adaptive filtering algorithms such 
as Least Mean Square (LMS) filters, RLS is popular because 
of its fast convergence and low steady-state error. 

In many applications it is natural to constrain the predictor 
coefficients to be sparse. In such cases the adaptive FIR 
prediction filter is a sparse system: only a few of the impulse 
response coefficients are non-zero. Sparse systems can be 
divided into general sparse systems and group sparse systems 
|[T], Unlike a general sparse system, whose impulse re- 
sponse can have arbitrary sparse structure, a group sparse sys- 
tem has impulse response composed of a few distinct clusters 
of non-zero coefficients. Examples of group sparse systems 
include specular multipath acoustic and wireless channels [3], 
1 4) and compressive spectrum sensing of narrowband sources 

The exploitation of sparsity to improve prediction perfor- 
mance has attracted considerable interest. For general sparse 
systems, the ii norm has been recognized as an effective 
promotor of sparsity [61, 17|. In particular, £1 regularized LMS 
||2), ||8) and RLS ||5),||ig algorithms have been proposed for 
for sparsification of adaptive filters. For group sparse systems, 
mixed norms such as the £1^2 norm and the £1.00 norm have 
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been applied to promote sparsity in statistical regression pT|- 
p3| , commonly referred to as the group lasso, and sparse 
signal recovery in signal processing and communications 
| fT4| . However, most of the proposed estimation algorithms 
operate in the offline mode and are not designed for time 
varying systems and online prediction. This is the motivation 
of our work. 

In this paper, we propose a RLS method penalized by the 
^1.00 norm to promote group sparsity, called the recursive £1.00 
group lasso. Our recursive group lasso algorithm is suitable for 
online applications where data is acquired sequentially. The 
algorithm is based on the homotopy approach to solving the 
lasso problem and is an extension of p3)-|[T7| to group sparse 
systems. 

The paper is organized as follows. Section II formulates 
the problem. In Section III we develop the homotopy based 
algorithm to solve the recursive £1,00 group lasso in an online 
recursive manner. Section IV provides numerical simulation 
results and Section V summarizes our principal conclusions. 
The proofs of theorems and some details of the proposed 
algorithm are provided in Appendix. 

Notations: In the following, matrices and vectors are de- 
noted by boldface upper case letters and boldface lower case 
letters, respectively; {■)'^ denotes the transpose operator, and 
II • 111 and II • II 00 denote the £1 and £^0 norm of a vector, 
respectively; for a set A, \A\ denotes its cardinality and (f> 
denotes the empty set; x_4 denotes the sub-vector of x from 
the index set A and Hab denotes the sub-matrix of R formed 
from the row index set A and column index set B. 

II. Problem formulation 

A. Recursive Least Squares 

Let w be a p-dimensional coefficient vector. Let y be an 
n-dimensional vector comprised of observations {yjYj^i- Let 
{xj}"^j^ be a sequence of p-dimensional predictor variables. 
In standard adaptive filtering terminology, yj, Xj and w are 
the primary signal, the reference signal, and the adaptive filter 
weights. The RLS algorithm solves the following quadratic 
minimization problem recursively over time n = p,p + 1, . . .: 

n 

w„ = argmin^7"~-'(?/j - w^Xj)^, (1) 
^ i=i 

where 7 e (0, 1] is the forgetting factor controlling the trade- 
off between transient and steady-state behaviors. 

To serve as a template for the sparse RLS extensions 
described below we briefly review the RLS update algorithm. 
Define R„ and r„ as 

n 

Rn = $]7""^x,xJ (2) 
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Fig. 1. Examples of (a) a general sparse system and (b) a group-sparse 
system. 



and 



The solution w„ to ([T]i can be then expressed as 

The matrix R„ and the vector r„ are updated as 

R„ = 7R„-i +x„x^, 

and 

Applying the Sherman-Morrison- Woodbury formula |[lF 

.-1 -1t>-1 



(4) 



where 



and 



gn -E^Ti— l-'^" 
1 



(5) 
(6) 

(7) 



Substituting Q into Q, we obtain the weight update 1 19| 

w„ = w„_i + a.„g„e„, (8) 

where 



(9) 



Equations (|5]l-(|9]) define the RLS algorithm which has com- 
putational complexity of order 0{p^). 



B. Non-recursive ii^oo group lasso 

The £i oo group lasso is a regularized least squares approach 
which uses the £i oo niixed norm to promote group- wise sparse 
pattern on the predictor coefficient vector The ii^oo norm of 
a vector w is defined as 



M 



W 1, 



where {Gm}m=i ^ group partition of the index set Q = 
{l,...,p}, i.e., 

M 

IJ Gm = Q, Gm n Gm' = (/i if TO ^ m', 

m— 1 

and wg^ is a sub-vector of w indexed by Gm - The oo norm 
is a mixed norm: it encourages correlation among coefficients 
inside each group via the £oo norm within each group and 
promotes sparsity across each group using the £i norm. The 
mixed norm ||w||i^oc is convex in w and reduces to ||w||i 
when each group contains only one coefficient, i.e., \Gi\ = 

\G2\ = ---^\Gm\^1. 

The ^1 oo group lasso solves the following penalized least 
squares problem: 



1 

arg mill - 



w^Xj)2 + All will, 



(10) 



(3) where A is a regularization parameter. Eq. ( 10 1 is a convex 



problem and can be solved by standard convex optimizers 



or path tracing algorithms |12|. Dkect solution of (lOi has 
computational complexity of 0{p'^). 

C. Recursive ii^ao group lasso 

In this subsection we obtain a recursive solution for ( fTO] ) 
that gives an update w„ from w„_i. The approach taken is 
a group-wise generalization of recent works p3) , |[T6l| that 
uses the homotopy approach to sequentially solve the lasso 
problem. Using the definitions Q and (j3]l, the problem (10 1 
is equivalent to 



w„ = argmin ^w^R„w - w^r„ + A||w||i,o 

w Z 



arg min (7R-n-i + x^x„ I w 



(11) 



w {jr 

) + ^l|w| 



Let f{/3,X) be the solution to the following parameterized 
problem 

■gmin^w'^ (7R„_i +;3x„x^) w 

vr Z (12) 

- w^(7r„_i + /3x„y„) + A||w||i,oo 
where /3 is a constant between and 1. w„ and w„_i of 



/(/3,A) 



problem (111 can be expressed as 
and 

w„-/(l,A). 

Our proposed method computes w„ from w„_i in the follow- 
ing two steps: 

Step 1. Fix /3 = and calculate /(0,A) from /(0,7A). This 
is accomplished by computing the regularization path between 
7A and A using homotopy methods introduced for the non- 
recursive £1 00 group lasso. The solution path is piecewise 
linear and the algorithm is described in |12|. 
Step 2. Fix A and calculate the solution path between /(O, A) 
and /(I. A). This is the key problem addressed in this paper. 
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Fig. 2. Illustration of the partitioning of a 20 element coefficient vector w 
into 5 groups of 4 indices. The sets V and Q contain the active groups and 
the inactive groups, respectively. Within each of the two active groups the 
maximal coefficients are denoted by the dark red color. 



To ease the notations we denote x„ and y„ by x and y, 
respectively, and define the following variables: 



R(/3) =7R„_i+/3xx^ 
r(/3) = 7r„_i + /^xy. 
Problem ^V2\ is then 



(13) 
(14) 



/(/3,A) = argniin-w^R(/3)w-w^r(/3) + A||w||i,oo. (15) 

w Z 



In Section III we will show how to propagate /(O, A) to /(I, A) 



using the homotopy approach applied to ( [T5] l. 

III. Online homotopy update 

A. Set notation 

We begin by introducing a series of set definitions. Figure 
|2] provides an example. We divide the entire group index set 
into V and Q, respectively, where V contains active groups 
and Q is its complement. For each active group m £ V, we 
partition the group into two parts: the maximal values, with 
indices Am, and the rest of the values, with indices Bm- 



and 



Am — arg max | | , to S V, 



Bm Qm Am • 



The set A and B are defined as the union of the Am and Bm 
sets, respectively: 



A^\J Am, S = U S„ 



niGV 



Finally, we define 



and 



C= D Qm- 

Cm — Qin ri C 



B. Optimality condition 



The objective function in ( 15 i is convex but non-smooth 



as the ii^oo norm is non-differentiable. Therefore, problem 



( 15 1 reaches its global minimum at w if and only if the sub- 



differential of the objective function contains the zero vector 
Let 9||w||i,oo denote the sub-differential of the oo norm at 
w. A vector z G 9||w||i only if z satisfies the following 
conditions |jl2), |[I4): 



llz.A„Jli = l,meP, 
sgn(z^„) = sgn(w^„ 
zg = 0, 

||zc„||i < i,™e Q, 



,meV, 



(16) 
(17) 
(18) 
(19) 



where A^B^C^V and Q are /3-dependent sets defined on w 
as defined in Section BlI-AI 

For notational convenience we drop /3 in R(/3) and r(/?) 
leaving the /^-dependency implicit. The optimality condition 
is then written as 



Rw — r + Az = 0, z G 9||w||i., 



(20) 



As wc = and zg = 0, (20 1 implies the three conditions 

R.^.AW^ + R^ewe - + Az^ = 0, (21) 
'^BA^A + Reewe - = 0, (22) 
'RcA^A + Rcewe - rc + Azc = 0. (23) 

The vector w_4 lies in a low dimensional subspace. Indeed, 
by definition of Am, if \Am\ > 1 

\wi\ = \wi>\, G Am- 

Therefore, for any active group m E V, 

= SA,„am (24) 

where 

= ||wg,„||oo, 

and 

sa = sgn (w^) . 
Using matrix notation, we represent (j24]l as 

w_A = Sa. (25) 

where 



S = 



(26) 



is a X IT-" I sign matrix and the vector a is comprised of 

am, m eV. 

The solution to ([TSj can be determined in closed form if the 
sign matrix S and sets {AjBjCjV, Q) are available. Indeed, 
from ^ and ([l7j 

S^z^ = 1, (27) 



where 1 is a jT'l x 1 vector comprised of I's. With (25 1 and 
( [ZT] ), ( [2T| i and ( [22] l are equivalent to 



S^R^^Sa + S'^R^ewg - S^r^ + Al = 0, 
Rg^Sa + RgBWe - rg = 0. 



(28) 
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Therefore, by defining the (a.s. invertible) matrix 



H 



R 



and 



BB 

a 



(29) 



(30) 



(28 I is equivalent to Hv = b — Ae, where e — (1 , ) , so 



that 



V = H-i(b- Ae). 



(31) 



As w_4 = Sa, the solution vector w can be directly obtained 



from V via (30i. For the sub-gradient vector, it can be shown 
that 

(32) 
(33) 

(34) 



and 



ze 

Azc = rc - (RcA^ Rcb) v. 



C. Online update 

Now we consider ( 15 1 using the results in III-B Let /3o and 
Pi be two constants such that /3i > /Jq- For a given value of 
/3 e [^o,/3i] define the class of sets S ^ {A,B,C,'P,Q) and 
make /3 explicit by writing S{P). Recall that S{(5) is specified 
by the solution /(/3, A) defined in (19). Assume that 5(/3) does 
not change for /? g [/3o, /3i]. The following theorem propagates 
/(/3o, A) to /(/3i, A) via a simple algebraic relation. 

Theorem 1. Let /3o and /3i be two constants such that /3i > /?□ 



and for any j3 € [/3o,/3i] the solutions to (751 share the same 
sets S — {A, B, C, V, Q). Let v' and v be vectors defined as 
A) and f{/3o,X), respectively. Then 



1 + ^jiPi 



(35) 



and the corresponding sub-gradient vector has the explicit 
update 



Az^ = A 



and 



Az^ = Azc 



1 + alPi 



1 + 



{y - y) {^A - (R-aaS R^e)g} 

(36) 



(y - y) {xc - (Rc^S RcbM 



(37) 

where R = R(0) as defined in [13), (x,2/) is the new sample 



as defined in \13\ and (14\, the sign matrix S is obtained from 



the solution at (3 = (3q, Hq is calculated from (29) using S 
and R(0), and d, u, y and are defined by 



= d^g. 



(38) 

(39) 
(40) 
(41) 



The proof of Theorem [T] is provided in Appendix A. 
Theorem [T] provides the closed form update for the solution 
path f{f3Q,X) /(/3i,A), under the assumption that the 
associated sets S{f3) remain unaltered over the path. 

Next, we partition the range /3 G [0, 1] into contiguous 
segments over which S(/3) is piecewise constant. Within each 
segment we can use Theorem 1 to propagate the solution from 
left endpoint to right endpoint. Below we specify an algorithm 
for finding the endpoints of each of these segments. 

Fix an endpoint /3o of one of these segments. We seek a 
critical point (5i that is defined as the maximum f3i ensuring 
S{(3) remains unchanged within [/3o,/3i]. By increasing /3i 
from /3o, the sets S{/3) will not change until at least one of 
the following conditions are met: 
Condition 1. There exists i G A such that z'^ = 0; 
Condition 2. There exists i £ B,n such that = a^; 
Condition 3. There exists m E V such that a'„ = 0; 
Condition 4. There exists m e Q such that ||z^^||i — 1. 
Condition 1 is from ( [T7] i and ( fTS] !, Condition 2 and 3 are based 
on definitions of A and V, respectively, and Condition 4 comes 
from (16 1 and (19i. Following |[T2|, pO), the four conditions 



can be assumed to be mutually exclusive. The actions with 
respect to Conditions 1-4 are given by 
Action 1. Move the entry i from A to B: 

A^ A- {i},B ^ B\J{i}] 
Action 2. Move the entry i from B to A: 

A^ A\J{i},B^B-{i}] 
Action 3. Remove group m from the active group list 

V - {m}, Q ^ Q U {m}, 
and update the related sets 

Action 4. Select group m 

V ^VU {to}, Q- {m}, 
and update the related sets 

By Theorem [T] the solution update from to /3i is in 
closed form. The critical point of (3i can be determined in a 
straightforward manner (details are provided in Appendix B). 
Let fi^i\k = 1,...,4 be the minimum value that is greater 
than /3o and meets Condition 1-4, respectively. The critical 
point /3i is then 



{k) 



Bi = min B\' 

fc=l....,4 



D. Homotopy algorithm implementation 

We now have all the ingredients for the homotopy update 
algorithm and the pseudo code is given in Algorithm [T] 

Next we analyze the computational cost of Algorithm[T] The 
complexity to compute each critical point is summarized in 
Table|I] where N is the dimension of Hq. As iV = [T'I + |S| < 
1^1 + is upper bounded by the number of non-zeros in 



5 



Algorithm 1: Homotopy update from /(O, A) to /(I, A). 

Input : /(0,A),R(0),x,y 
output: /(I, A) 

Initialize /3o = 0, = 0, R = R(0); 

Calculate {A,B,C^'P iQ) and (v,Az^,Azc) from 

/(0,A); 

while /3o < 1 do 

Calculate the environmental variables 
(S,Ho,d,g,2/,cr|^) from f{l3o,X) and R; 
Calculate {f3[''^}t^i that meets Condition 1-4, 
respectively; 

Calculate the critical point f3i that meets Condition 

k^,: k^, — arg miiife /?^'^'' and jSi — P^''^; 
if /3i < 1 then 



Update (v, Az_4,Azc) using (35 i, (36i and (37i; 
Update {A, B, C, V, Q) by Action 

/3o = 

else 
I break; 
end 
end 

/3i = i; 

Update (v, Az^,Azc) using (35 i; 
Calculate /(I, A) from v. 



the solution vector The vector g can be computed in 0{N'^) 



time using the matrix-inverse lemma 1 1 8 1 and the fact that, 
for each action, Hq is at most perturbed by a rank-two matrix. 
This implies that the computation complexity per critical point 
is 0{pma,x{N, logp}) and the total complexity of the online 
update is 0(^2 ■ pmax{N,logp}), where fc2 is the number 
of critical points of /3 in the solution path /(O, A) — > /(I, A). 
This is the computational cost required for Step 2 in Section 

ED 

A similar analysis can be performed for the complexity of 
Step 1, which requires 0{ki • pmax{A^, logp}) where fci is 
the number of critical points in the solution path /(0,7A) — > 
/(O, A). Therefore, the overall computation complexity of the 
recursive £i,oo group lasso is 0{k ■ pmax{N,\ogp}), where 
k = ki + k2, i.e., the total number of critical points in the 
solution path /(0, 7A) ^ /(O, A) ^ /(I, A). 

An instructive benchmark is to directly solve the n-samples 



problem (12i from the solution path /(l,cx)) (i.e., a zero 
vector) /(I, A) | [T2) , without using the previous solution 
w„_i. This algorithm, called iCap in fT2'|, requires 0{k' ■ 
pma,x{N, \ogp}), where k' is the number of critical points in 
/(I, 00) — > /(I, A). Empirical comparisons between k and k', 
provided in the following section, indicate that iCap requires 
significantly more computation than our proposed Algorithm 

m 



IV. Numerical simulations 

In this section we demonstrate our proposed recursive ii^ao 
group lasso algorithm by numerical simulation. We simulated 



g = H-M 






0(\A\N) 




n(\r \ N^ 




Oi\A\) 




0(|B|) 




Oi\V\) 




0(|C|log|C|) 



TABLE I 

Computation costs of online homotopy update for each 
critical point. 
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Fig. 3. Responses of the time varying system, (a): Initial response, (b): 
Response after the 200th iteration. The groups for Algorithm 1 were chosen 
as 20 equal size contiguous groups of coefficients partitioning the range 
1,...,100. 



zero mean Gaussian noise and w* is a sparse p = 100 element 
vector containing only 14 non-zero coefficients clustered be- 
tween indices 29 and 42. See Fig. [3] (a). After 200 time units, 
the locations of the non-zero coefficients of is shifted to 
the right, as indicated in Fig. [3] (b). 

The input vectors were generated as independent identically 
distributed Gaussian random vectors with zero mean and iden- 
tity covariance matrix, and the variance of observation noise 
Vj is 0.01. We created the groups in the recursive oo group 
lasso as follows. We divide the 100 RLS filter coefficients 
w into 20 groups with group boundaries 1, 5, 10, . . ., where 
each group contains 5 coefficients. The forgetting factor 7 
and the regularization parameter A were set to 0.9 and 0.1, 
respectively. We repeated the simulation 100 times and the 
averaged mean squared errors of the RLS, sparse RLS and 
proposed RLS shown in Fig. |4] We implemented the standard 
RLS and sparse RLS using the £1 regularization, where the 
forgetting factors are also set to 0.9. We implemented sparse 
RLS 1 15 1 by choosing the regularization parameter A which 
achieves the lowest steady-state error, resulting in A — 0.05. 

It can be seen from Fig. [4] that our proposed sparse RLS 
method outperforms standard RLS and sparse RLS in both 
convergence rate and steady-state MSE. This demonstrates 
the power of our group sparsity penalty. At the change point 
of 200 iterations, both the proposed method and sparse RLS 
of fl5j show superior tracking performances as compared to 
the standard RLS. We also observe that the proposed method 
achieves even smaller MSE after the change point occurs. This 
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Fig. 4. Averaged MSB of the proposed algorithm, RLS and recursive lasso. 
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Fig. 5. Averaged number of critical points for the proposed recursive method 
of implementing ii^oo lasso and the iCap (12| non-recursive method of 
implementation. 



is due to the fact that the active cluster spans across group 
boundaries in the initial system (Fig. [3] (a)), while the active 
clusters in the shifted system overlap with fewer groups. 

Fig. [5] shows the average number of critical points (ac- 
counting for both trajectories in /3 and A) of the proposed 



algorithm, i.e., the number k as defined in Section III-D 
As a comparison, we implement the iCap method of ||12), 
a homotopy based algorithm that traces the solution path only 
over A. The average number of critical points for iCap is 
plotted in Fig.js] which is the number k' in Section III-D Both 
the proposed algorithm and iCap yield the same solution but 
have different computational complexities proportional to k 
and k', respectively. It can be seen that the proposed algorithm 
saves as much as 75% of the computation costs for equivalent 
performance. 

V. Conclusion 

In this paper we proposed a £1^00 regularized RLS algorithm 
for online sparse linear prediction. We developed a homotopy 
based method to sequentially update the solution vector as new 
measurements are acquired. Our proposed algorithm uses the 
previous estimate as a "warm-start", from which we compute 
the homotopy update to the current solution. The proposed 
algorithm can process streaming measurements with time 



varying predictors and is computationally efficient compared 
to non-recursive group lasso solvers. Numerical simulations 
demonstrated that the proposed method outperformed the 
standard and £1 regularized RLS for identifying an unknown 
group sparse system, in terms of both tracking and steady-state 
mean squared error. 

The work presented here assumed non-overlapping group 
partitions. In the future, we will investigate overlapping groups 
and other flexible partitions [ |21J . 

VI. Appendix 
A. Proof of Theorem |7] 

We begin by deriving ( [35] l. According to (31 1, 

v' = H'-i(b'-Ae'). (42) 

As S and {A,B,C,V^ Q) remain constant within [/3o,/?i], 

e' = e, (43) 

b' = b + 5dy, (44) 

and 



where 



H' = H + 5Ad^ 



5 = Pi- (3o, 



H and b are calculated using S within [/3o, /3i] and R(/3o) and 
r(^o)^ respectively. We emphasize that H is based on R(/3) 
and is different from Hq defined in Theorem [T] According to 
the Sherman-Morrison-Woodbury formula, 

S 



H' 



H 



1 + a'^S 



(H-M)(H-M)^ 



(45) 



where = d^H~M. Substituting (43 1, (44 1 and (45 1 into 



( 42 1, after simplification we obtain 

5 



H 



1 + 0-25 



{H-^d){H-^df (b + Sdy - Ae) 



1 + 0-2,5 
S 



S 



H Md'H^(b Ae) 



^2^2 



U^dy 



1 + a'^o 

(46) 

where y = d^v as defined in ( 40 1. 

Note that H is defined in terms of R(/3o) rather than R(0) 
and 

H = Ho + /3ocid^, 



so that 



H 



Hn 



/?0 



(47) 



where g and ajj are defined by (39i and (41 1, respectively. 



As ajj 



H M = Hr^M 



(48) 
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Accordingly, 



H 



1 + (^IPG 



•h 



'H 



1 + (ThPo ' 



(49) 



Substituting (48 i and (p9]l to d46]), we finally obtain 



-(2/ - y)g = V + f] fl [y - j})g. 



1 + ajjlii "'^ 1 + alPi ' 

Equations ([36] l and ( [37] i can be established by direct sub- 
stitutions of p5| l into their definitions ( |32| l and ( (34] l and thus 
the proof of Theorem [T] is complete. 

B. Computation of critical points 

For ease of notation we work with p, defined by 

P - T , ^2 o ■ (50) 

It is easy to see that over the range /3i > /3o, p is monotonically 
increasmg m (0, l/cr|j). Therefore, JsOb can be inverted by 



Pi 



1 ? ' 



(51) 



where p e (0, l/cr|j) to ensure /3i > /3o- 

Suppose we have obtained p^^\k = 1,...,4, /Sj;'^-' can be 



calculated using (51 1 and the critical point /3i is then 



j3i = min /Sl'^''. 

fc=l,...,4 

We now calculate the critical value of p for each condition 
one by one. 

1) Critical point for Condition 1: Define the temporary 
vector 

t.A = (y - {X^ - (R^.aS R^i3)g}- 

According to ( [36| , 



Condition 1 is met for any p 



.(1) 



such that 



17 



,i <^ A. 



Therefore, the critical value of p that satisfies Condition 1 is 



p'^^^ = mm • 



{pr'|zeApr'e(0,i/a|,)} 



2) Critical point for Condition 2: By the definition pUf , v 
is a concatenation of a„i and we,„ ,m £ V: 



= (^(a„i)™e-p,we^,...,we|^J , (52) 
where {am)me'P denotes the vector comprised of a„i, m ^V. 



Now we partition the vector g in the same manner as ( 52 1 and 
denote t„j and u„j as the counter part of q;,„ and wg^ in g, 
i.e., 

T ( I \ T T 

g — UT„jj,„g-p,Ui ,...,U|p| 



Eq. ( 35 I is then equivalent to 



(53) 



and 



+ pUm,i, 



where Um.i is the z-th element of the vector u,„. Condition 2 
indicates that 



qL — ±w'r. 



and is satisfied if p 



rrn.i 



PmA or p 



Pm.i^' where 



+ T„ 



Therefore, the critical value of p for Condition 2 is 



(2) • r (2±) 

P ' = mm|p;„^^ 



TO G T', i = 1. 



I'-'nil'Pm.i ^ 



(0,1/^1,)} 



yields p ~ p\ ' determined by 



3} Critical point for Condition 3: According to (53 1, q;J„ 
errr 

(3) 



, TO, S "P, 



and the critical value for p^''^ is 



,TO e 75,^^^) e 



4j Critical point for Condition 4: Define 

tc = (y - {xc - (RcaS Rce)g}. 



Eq. (37 1 is then 

Azg^^ = Azc,„ + ptc„ 
and Condition 4 is equivalent to 

+ Azil = A 



(54) 



iec„ 



To solve ( [54] i we develop a fast method that requires complex- 
ity of 0{N \ogN), where N — \Cm\- The algorithm is given 
in Appendix C. For each to e Q, let p^m be the minimum 



positive solution to (54i. The critical value of p for Condition 
4 is then 



P 



(4) 



min {pW|TOeQ,p(t^e (0,1/^1,)} 



C. Fast algorithm for critical condition 4 

Here we develop an algorithm to solve problem ( [54| . 
Consider solving the more general problem: 



N 

E 



ai\x - Xi\ = y, 



(55) 



where and Xi are constants and > 0. Please note that 
the notations here have no connections to those in previous 
sections. Define the following function 



h{x) 



N 

E 



Qi \x Xi I . 



The problem is then equivalent to finding h^^{y), if it exists. 

An illustration of the function h{x) is shown in Fig. |6] 
where ki denotes the slope of the ith segment. It can be shown 
that h{x) is piecewise linear and convex in x. Therefore, the 
equation ( [55] l generally has two solutions if they exist, denoted 
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Fig. 6. An illustration of the fast algorithm for critical condition 4. 



as a;min and Smax- Based on piecewise linearity we propose a 



search algorithm to solve (55 1. The pseudo code is shown in 



Algorithm |2] and its computation complexity is dominated by 
the sorting operation which requires 0{N\ogN). 



Algorithm 2: Solve x from X]il=i "^'N 



Input : {a,,a;J^i, y 

output. Xjnin 1 ^max 

Sort {xi}fLi in the ascending order: xi < X2 < ... < x^; 
Re-order {ai}fL^ such that corresponds to x^; 

Set ko = -J2i=i(^i' 
for « = 1, N do 

end 

Calculate hi = X]i!=2 '^iNi ~ 
for i = 2,...,N do 

I /i, = hi^i + k^^iixi - 
end 

if mini ki > y then 
No solution; 
Exit; 
else 

it y > hi then 

I a;,ni„ = + (y - /ii)/fco; 
else 

Seek j such that y e [hj, hj^i]; 

2^min = + (j/ ~ 

end 

it y > hff then 

I a^max = XJV + (y - hN)/kN; 

else 

Seek j such that i/ e [/ij-i, /ij]; 
a^max = Xj^i + (y - hj_i)/kj_i; 

end 
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